Monte Carlo Study of the Axial Next-Nearest-Neighbor Ising Model 
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The equilibrium phase behavior of microphase-forming systems is notoriously difficult to obtain 
because of the extended metastability of the modulated phases. We develop a simulation method 
based on thermodynamic integration that surmounts this problem and with which we describe 
the modulated regime of the canonical three-dimensional axial next-nearest-neighbor Ising model. 
Equilibrium order parameters are obtained and the critical behavior beyond the Lifshitz point is 
examined. The absence of widely extended bulging modulated phases illustrates the limitations of 
various approximation schemes used to analyze microphase-forming models. 

PACS numbers: 64.60.Cn, 64.60.F-,05.10.Ln,75.10.-b 
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Microphases self-assemble in systems with competing 
short-range attractive and long-range repulsive interac- 
tions, irrespective of the physical and chemical nature 
of these interactions Microphases are the frustrated 
equivalent of gas-liquid coexistence for purely attract- 
ing particles. Periodic lamellae, cylinders, clusters, etc. 
are thus observed in a variety of systems, such as multi- 
block copolymers Q, oil- water surfactant mixtures Q, 
charged colloidal suspensions Q, and magnetic materi- 
als [5|. Although the modulated organization is spon- 
taneous, obtaining detailed morphological control is no- 



external fields [7|, or 
are usually necessary 



toriously difficult. Annealing [6 
complex chemical environments 
to order diblock copolymers. Mesoscale periodic textures 
have found some technological success as thermoplastic 
elastomers Q and nanostructure templates Q, but un- 
derstanding how to tune and stabilize microphases is es- 
sential to broadening their material relevance. 

Because experimental systems provide only limited mi- 
croscopic i nsig ht into microphase formation, a number of 
lattice [10l4l3l | and free-space EMI models have been 
put forward. Grasping the equilibrium properties of these 
models is necessary to resolve the problems surround- 
ing the non-equilibrium assembly of microphases 17 - 19| ■ 
Though the modulated regime is a central feature of these 
systems, microphases have not been accurately charac- 
terized in any of them. Even for simple models, ap- 
proximate theoretical frameworks offer only limited as- 
sistance, and treating microphases with computer sim- 
ulations is so far an unresolved problem |20l 1 2 1 1 . In 
this Letter, we overcome this last issue by developing 
a free-energy integration method for modulated phases. 
We use this method to determine the phase diagram 
of the microphase-forming three-dimensional (3D) axial 
next-nearest-neighbor Ising (ANNNI) model, which has 
reached textbook status [HI [22|, but whose character- 
istic modulated behavior is still not completely under- 
stood. The resulting phase information allows us to as- 
sess the validity of competing approximate treatments 
and to better understand the phenomenology of related 
experimental systems. 



The ANNNI model was introduced nearly half a cen- 
tury ago to explain "helical" magnetic order in heavy 
rare-earth metals 0,1 [23M25I] . Its Hamiltonian on a sim- 
ple cubic lattice for spin variables s$ = ±1 
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favors alignment of nearest- neighbor pairs but frus- 

trates long-range order with relative strength k > for z- 
axial next-nearest-neighbor pairs The coupling con- 
stant J determines the temperature T scale with Boltz- 
mann's constant fcs set to unity for convenience. The 
ANNNI model can only be solved exactly in one dimen- 
sion [26[ , but some of its higher-dimensional features are 
nonetheless well understood. In 3D, the topography of 
the T-k phase diagram involves three regions that join 



together at a multicritical Lifshitz point [27|: at high T 
the system is paramagnetic; at low T and k it is ferro- 
magnetic; at low T and for sufficiently high k modulated 
layered phases form (25|. The ANNNI paramagnetic- 
modulated (PM) transition beyond the Lifshitz point is 
thought to be part of the XY universality class [28[ . For 
k < 1/2 the T = ground state is ferromagnetic, and for 
k > 1/2 it is the layered antiphase ("two- up- two-down"). 
The sequence of commensurate phases springing from the 
multiphase point at T = and k = 1/2, the structure 
combination branching processes at low T, and the pos- 
sible occurrence of incommensurate phases are also note- 
worthy features of the model [29| . 

In order to detail the phase behavior, approximate the- 
oretical treatments, including high- and low-temperature 
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simulations 




mean-field 
been used. 



32j, |33(, and 



have been used. Monte Carlo 
have reliably determined the 
paramagnetic-ferromagnetic transition up to the Lifshitz 
point [40(, but accurately locating transitions to and 
within the modulated regime has remained elusive. Even 
within the subset of periodic phases commensurate with 
the finite lattice, high free-energy barriers need to be 
crossed on going from one modulated phase to another. 
Patterns with a metastable nearby periodicity thus per- 
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FIG. 1: (Color online) Simulation results for k = 0.7. (Top) 
Energy and free energy per spin for modulations ranging from 
q = 1/4 (antiphase) to q c = 0.1917 at melting. The PM tran- 
sition T c = 3.988(1) (dashed line) is extracted from suscep- 
tibility measurements (Fig. [2}. (Bottom) Equilibrium devil's 
staircase and generalized magnetization m(q). The power- law 
decay of m(q) with /3 = 0.34(4) is superimposed (line). (In- 
set) Snapshot of the antiphase with differently shaded beads 
for Si — ±1. 



sist for very long times 24, 26, 38]. Traditional sim- 
ulation methodologies that facilitate ergodic sampling 
of phase space by passing over such barriers, notably 
parallel tempering and cluster moves, are of limited 
help in microphase-forming systems. Because the equi- 
librium periodicity varies with temperature, sampling 
higher temperatures leaves the system in a modulated 
phase with the wrong periodicity; because of the high 
free-energy barriers between modulated phases and the 
lack of simple structural rearrangements for sampling dif- 
ferent modulations, the efficiency of cluster moves is lim- 
ited. 

We develop a simulation method based on free-energy 
integration to treat microphases. The free energy of mod- 
ulated phases allows us to compare the stability of differ- 
ent periodic patterns and to reliably capture phase tran- 
sitions. Some aspects of the procedure are part of the 
standard numerical toolkit 41| , but additional specifica- 
tions are in order. For a given k, T, and wave num- 
ber modulation q, we first calculate the absolute free 
energy F of g-modulated lamellae at a nearby refer- 
ence temperature To, and then thermally integrate the 
energy per spin E/N from To to T. In the spirit of 



Refs. 42J, |43j , the Kirkwood integration begins from de- 
coupled spins under an oscillatory sinusoidal field with 
Hamiltonian Ho = —E>o X)t=i s « sm(2nqzi + 0o), where a 
small phase angle (f>o is added to prevent the lattice sites 
from overlapping with the zeros of the field. A scaling 
field Bo sufficiently strong to avoid melting is necessary 
for the reversibility of the integration scheme. The high 
free-energy barriers between the neighboring commensu- 
rate periodic patterns would also make phase transitions 
highly unlikely even if sections of the path are formally 
metastable [26}]. Similarly a sinusoidal reference state is 
valid even if the layer profile squares at low T (26[, be- 
cause there is no phase transition along the integration 
path. We perform constant T Monte Carlo (MC) simula- 
tions on a periodic lattice with N = L x L y L z = 40 2 x 240, 
unless otherwise noted. Wave numbers q — n/L z for in- 
teger values of n keep modulations commensurate with 
the lattice, which leaves open the problem of incommen- 
surate phases. Phase-space sampling gains in efficiency 
when single-spin flips are complemented with MC moves 
that take advantage of phase symmetries. In the modu- 
lated phases, layer exchanges allow for thickness fluctu- 
ations and lattice drifts sample the external field; in the 
paramagnetic phase, cluster moves accelerate sampling 
in the critical region [ijj. For To reference integrations, 
up to 10 5 MC moves (N attempted flips) are performed 
after 5 x 10 4 MC moves of preliminary equilibration. 

The smooth and extended energy curves of the dif- 
ferent modulations are characteristic of the long-lived 
metastable nature of the periodic phases (Fig. [T]). Even 
over relatively long simulation times, metastable systems 
do not relax to their equilibrium periodicity. Thorough 
sampling is possible without any modulation change if 
the L x L y cross-section is sufficiently large. The energy 
gap between neighboring phases for g's commensurate 
with the simulation box reflects the limited choice of 
modulations on a finite lattice. In an infinite periodic 
system, where all rational modulations are valid but ir- 
rational q's are excluded, the gap becomes infinitely small 
because rational numbers are dense on the real axis [38| . 
Although they appear to join together smoothly on the 
scale of Fig. [TJ neighboring free-energy curves intercept. 
The intercept identifies the transition temperature be- 
tween two modulated phases with an accuracy that vastly 
surpasses previous simulation approaches 
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ure Q] illustrates the profile of the devil's staircase for 
k = 0.7 [34{. The rate of change of the equilibrium q 
accelerates upon cooling. The predicted discontinuity of 
the function before reaching the antiphase should make 
the staircase "harmless" , but the current numerical ac- 
curacy is insufficient to distinguish this scenario from the 
"devil's last step" [13, H]- 

The PM critical transition temperature T c , which is an- 
alytically well characterized El 51, is used to validate 
the simulation results. Because the heat capacity per spin 
C/N is at best only weakly divergent at T c (Fig. [5]), we 
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FIG. 2: Finite-size scaling of x(?c) at k = 0.7 around the PM 
transition with v = 0.60(3) and 7/1/ = 2.13(3). (a) Same for 
C/N with a/u = 0.18(2). (b) Simulation and series expansion 
7 compared with Ising (k < kl) and XY (k > k_l) exponents 
(dashed lines); kl = 0.270 result from Ref. [40(]. Rushbrookc 
and hyperscaling equalities are obeyed within error bars. 



also consider order parameters that are functions of the 
Fourier spin density s q = Yli=i Sie l2 ' nqZi and thus natu- 
rally capture modulations. In the paramagnetic phase, 
the z-axis static structure factor 5(g) = (s q S- q )/N grows 
upon cooling and di verg es at the critical wave number 
q c obtained in Fig. Q] j2(| [30[ • But the system-size diver- 
gence of 5(g) on the modulated side makes it ill-suited for 
determining T c in simulations. The generalized magneti- 
zation per spin m(q) = N^ 1 ^/ (s q )(§- q ) also causes prob- 
lems, because it averages to zero as the lattice drifts fl3j . 
To correct for this problem, we maximize the real com- 
ponent of s q with a phase shift for each configuration, 
before taking the thermal average. The resulting func- 
tion shows the characteristic power law m(q) ~ \T — T C \P 
decay (Fig. [T]). The transition is, however, most clearly 
identified from the generalized Binder cumulant [2l| (not 
shown) and the generalized susceptibility 

NT X {q) = (s q ~S- q )-{s q ){s- q ) ^NS{q)-N 2 m 2 {q) 1 (2) 

which diverges with system size x(g c ) ~ \T — T c |~ 7 
(Fig. [2|) as does x(0) at an Ising-like transition. The 
T c results are in very good agreement with the series ex- 
pansion 3(1 45 1, indicating that T c can be identified to 
within a part in a hundred using C from the standard 
lattice size. The resulting determination of the PM tran- 
sition (Fig. Ep i s also more reliable than the rare earlier 
MC results |26l |4o| . because of the larger system sizes 
used. 

We also examine the suggested XY character of the PM 
transition (H- The derivative of hx{S{q)/N) with J/T 
gives the correlation length divergence exponent is, while 
the exponent ratios a/u and 7/V are determined by finite- 
size scaling of C/N and x(g), respectively (Fig. [5]) [2l| . 




FIG. 3: (Color online) Lifshitz point (•) [401 an d simulation 
phase boundaries from x(lc) (□), C (x), and F (0 and dotted 
line) . High- [30l . 1451 ] and low-temperature (up to third order 
around the multiphase point) I3lj series expansions as well as 
effective-field [35( and TPVA [3a | results are indicated. Sta- 
bility region of phases q — 1/6 (left inset) and q = 1/5 (right 
inset) are compared with different theoretical approaches. 



The exponent ratios above kl, though consistent with 
each other, may suggest that the PM transition has a 
universality that is not of XY type. In particular, a has 
a positive sign and 7/1/ is significantly different from the 
XY value for the ratio, which may explain the discrep- 
ancy of the series expansion 7 results for large k (see 
Fig. M and caption) [30, SB | . 

More significantly, the approximate treatments, which 
capture the external boundaries of the modulated regime 
reasonably well 35|, 36} , qualitatively disagree on the in- 
ternal structure of that regime. On the one hand, as 
with the mean-field treatment [32I [33[ and the soliton 
approximation [34[, the effective-field method fills the 



modulated interior with exceptionally stable "simple pe- 
riodic" [3] bulging phases, such as the "three- up-three- 
down" q = 1/6 phase and the q = 1/5 phase [35(. On 
the other hand, the tensor product variational approach 
(TPVA) predicts rather narrow widths for the commen- 
surate phases [36| . The simulation results bulge less than 
is suggested by the first scenario. The rate of change of 
q with k and T slows in certain parts of the modulated 
regime, but all of the phases commensurate with the pe- 
riodic box are stable in turn. The stability range of the 
different modulated phases is overall fairly small and no 
exceptional stability is observed for the simple periodic 
phases q = 1/6 and q — 1/5 (Fig. [3]), unlike for q = 1/4. 
The q = 1/6 phase does bulge, but increasing the system 
size, which allows for a more refined q selection, results 
in a shrinking stability range (Fig. [3]), in opposition to 
the q = 1/4 phase whose stability range is system-size 
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independent. For the q = 1/5 phase, the range of stabil- 
ity increases slightly with k in simulation, which is also 
due to the finiteness of the lattice. It is possible that 
the reduced range of stability of these phases compared 
to the mean-field predictions be related to the relatively 
low roughening transition (T r — 2.445 47 1) in the cor- 
responding Ising model compared to the temperatures 
studied here. Further study is needed to clarify this 
point. The absence of widely extended bulging phases 
suggests that the lack of qualitative agreement between 
observations in magnetic systems, such as CeSb [jlHU, 
and the mean-field stability ranges is to be expected. The 
commensurate phases observed are those that are kineti- 
cally accessible on experimental time scales [26| or whose 
stability is due to corrections beyond simple spin mod- 
els. Neither effect suggests a preferable agreement with 
mean-field predictions. 

In this Letter we have presented a methodology for 
simulating layered microphases, but modulated assem- 
blies can exhibit a variety of other symmetries, under 
the control of an external magnetic field or by tun- 
ing the chemical potential in the corresponding lattice 
gas model. Generalizing the approach to other order 
types will greatly benefit the study of more elaborate 
microphase-forming systems and pave the way for stud- 
ies of the non-equilibrium microphase assembly, where 
most of the materials challenges lie. Generalization to 
frustrated quantum systems is also conceivable as long 
as the sign problem can be surmounted |49j . 
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edge ORAU and Duke startup funding. 
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